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Abstract 

We study a restricted- height version of the one-dimensional Oslo sandpile with conserved 
density, using periodic boundary conditions. Each site has a limiting height which can be 
either two or three. When a site reaches its limiting height it becomes active and may 
topple, loosing two particles, which move randomly to nearest-neighbor sites. After a site 
topples it is randomly assigned a new limiting height. We study the model using mean- 
field theory and Monte Carlo simulation, focusing on the quasi-stationary state, in which 
the number of active sites fluctuates about a stationary value. Using finite-size scaling 
analysis, we determine the critical particle density and associated critical exponents. 
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I. INTRODUCTION 



Sandpile models are paradigmatic examples of self-organized criticality (SOC) 
, a control mechanism that forces a system with an absorbing-state phase 
transition to its critical point js-S], without explicit tuning of control parameters {g]. 
SOC in a slowly driven sandpile corresponds to an absorbing-state phase transition 
in a model with the same local dynamics, but a fixed number of particles jsl, 



so-called conserved sandpiles |lOl. Il2|-|l4l|. Absorbing-state phase transitions arise in 
the context of spatial stochastic models, and correspond to a transition between an 
active, fluctuating phase, and an absorbing one, which allows no escape 

Sandpile models with probabilistic toppling rules, typified by the Manna model 
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16| . are commonly designated as stochastic sandpiles; their study has been 



central to establishing the connection between SOC and absorbing-state phase tran- 
sitions. An important stochastic model is the Oslo model inspired by experi- 
mental studies on rice piles. In this work we study a conserved version of the Oslo 
model, characterizing its absorbing-state critical point. 

An inconvenient feature of many sandpile models is the absence of an upper 
bound on the number of particles that may occupy a given site, which complicates 
theoretical approaches such as n-site approximations or continuum descriptions. 
This motivated the study of restricted sandpiles 



. In the present work we impose a 



height restriction on the conserved Oslo model. Since the symmetries and conserved 
quantities of the restricted and unrestricted models are the same, one expects, on 
the basis of experience with critical phenomena both in and out of equilibrium, that 
the models belong to the same universality class, as is indeed borne out for conserved 
versions of the Manna model 
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2l|. The symmetries here are limited to spatial 



translation and inversion, while the conservation law is that of particle number. This 
universality class has come to be known as the conserved directed percolation (CDP) 
class. Recently it was suggested that the critical behavior of conserved stochastic 
sandpiles in fact belongs to the (nonconserved) directed percolation class 22|, but 
further studies are required to verify this assertion. 

The remainder of this paper is organized as follows: In Section II we describe the 
model and in Section III develop a one-site mean-field approximation. Our numerical 



results are reported in Section IV and in Section V we present our conclusions. 
II. MODEL 

We study a restricted sandpile model with a variable height limit, defined on a 
lattice of NsUe = L'^ sites, where d is the dimension of space. The configuration is 
specified by the particle numbers Zi [i = 1, . . . , A^^jfe) at each site. Each site has 
a critical height Ci G {2,3}, with equal probability, such that Zi < Ci. A site with 
Zi = Ci is said to be active. Any configuration devoid of active sites is absorbing, 
i.e., it admits no escape. The dynamics of the model proceeds via toppling of active 
sites; each active site has a rate of unity to topple. (By a "rate of unity" we mean 
that the time unit is chosen such that if there are currently A^^ active sites, then the 
time increment associated with the next toppling is 1/Na-) In a continuous-time 
(sequential) dynamics, each active site has the same probability of being the next 
to topple. When a site, say i, topples, two particles are transferred from i to sites j 
and j', nearest neighbors of i. The two sites are chosen at random, independently, 
from the set of nearest neighbors, and so are not necessarily distinct; we refer to 
this procedure as an independent toppling rule. Due to the height restriction, any 
particle transfer that would result in a target site j having Zj > cj is rejected. (This 
means that the configuration Zi = Ci, Vz is also absorbing, since no particles can 
be transferred. The particle densities of interest in this study, however, remain far 
below the density associated with this configuration.) When site i loses a particle 
or particles due to toppling, a new limiting height q is selected, equal to 2 or 3, 
each with probability 1/2. Thus the dynamics has three stochastic elements: (1) the 
choice of the next site to topple; (2) the choice of target sites j and j' for particle 
transfers; (3) the choice of the new limiting height after a site topples. 

In practice the next site to topple is selected at random from a list of currently 
active sites, which must naturally be updated following each toppling event. The 
time increment associated with each toppling (whether particles are transferred or 
not) is At = 1/Na, with A^^ the number of active sites immediately prior to the 
event. 
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III. MEAN-FIELD THEORY 



The primary aim of the mean-field analysis is to obtain a preliminary idea of 
the phase diagram and (assuming the latter possesses a phase transition), an order 
of magnitude estimate of the critical point. We consider the simplest mean- field 
approach, known as the one-site approximation. At this level of approximation, 
there are seven possible states for a given site, where i = 0, 1, j represents 
the occupation number and j = 2, 3 denotes the limiting height, with associated 
probabilities denoted by Pij. Taking into account the conditions of normalization 

j=2 i=0 

and of fixed density, 

EE^^.. = c, (2) 

j=2 i=0 

there are only five independent variables at this level. 

We begin the analysis by listing the possible transitions between states in 
Fig. [TJ Each transition (at a given site, called the central site in this discussion), 
requires a specific configuration at the central site and at one or both of its nearest 
neighbors, and a certain redistribution of particles from the toppling site. (The local 
configuration and the choice of target sites, j and j', in the particle redistribution are 
statistically independent events.) In the one-site approximation, joint probabilities 
involving two or more sites are factorized. Denoting a joint two-site probability by 
Pij\ki, the one-site approximation uses the replacement Pij\ki PijPki, and similarly 
for three-site probabilities. 

To illustrate how the rates associated with these transitions are calculated, we 
discuss some examples. Consider first the transition (0|2) (1|2). The initial 
configuration must be either (02|22) or (03|23), that is, the central site must be 
vacant, have Zc = 2, and have an active neighbor. When the latter topples, exactly 
one particle must migrate to the central site. On a hypercubic lattice in d dimensions, 
each site has 2d nearest neighbors. Since the probability of exactly one particle 
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FIG. 1: Transitions between states of a single site. Transitions marked "X" are impossible 
and those ones related to diagonal elements are irrelevant. 



jumping to the central site is 2{l/2d)[l — {l/2d)], the rate (per site) of transitions 
of the kind (0|2) ^ (1|2) is 

2d — 1 

2rf^^P02(P22 + P33), (3) 

where the factor 2d represents the number of nearest neighbors. 

Consider next the transition (2|2) — (0|2), which can occur via two mutually 
exclusive paths. In one, both particles liberated when the central site topples at- 
tempt to migrate to the same neighbor, an event having probability l/Ad"^. In order 
for both particles to actually migrate, the difference Zc — z ai the target site must 
be greater than one. Thus the initial configurations for which this transition may 
occur are (00|23), (02|32), and (12|32). The transition rate for this path is 

^2d^P22{P02 + P03 + Pl3), (4) 

where the factor 1/2 represents the probability that the limiting height retains the 
value of 2 following the toppling event. In the other path, the two particles migrate 
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to distinct neighbors of the central site. The configurations that allow this transi- 
tion to occur are (2^22^1222), (2^23^1223), and (3+23^1323), where 2^ and 3+ denote, 
respectively, sites with z <2 and 2 < 3. Thus the transition rate for this path is 

2d — 1 

^^P22(At2 + P3t3)'. (5) 

Evaluating the rates of the remaining transitions, we find the equations that govern 
the probabilities Pjj at this level of approximation. The equations for the Pij are 



dPmit) 
dt 

dt 

dP23{t) 

dt 

dPi2{t) 
dt 



dt 



dP22{t) 

dt 



dP33{t) 

dt 



— P22(-Po2 + -P03 + -P13) H 73 — -P22(-P2t2 + -P3t3)^ — -^02 



4d 
2d -I 



P22{P2^2 + Pz^3f-^^P^ 



03 



Ad 

x(P22 + i^33), 

^P22(P02 + P03 + ^3) + - . - o^o. 

X(P22 + P33) 

^ [Po3 + 2{2d - l)Pi3 - {Ad - 1)P23] (P22 + P33) + + 2(2c/ - 1: 



(6) 



(7) 



xiP22 + P33)]iPl2 + P23)P 



33, (8) 
^^P02 - ^^^12) (^22 + P33) + ^ [1 + 2(2rf - 1)(P22 + P33)] 

XP22(P12 + P23) + ^ [P02 + P03 + Pl3 + i2d - l)(P2t2 + P3t3)'] ^33, 

(9) 

^^P03 - ^^^13) (^22 + P33) + ^ [1 + 2(2d - 1)(P22 + P33)] 



X 



P22(Pl2 + P23) + ^ [^02 + P03 + Pl3 + {2d - l)(P2t2 + P3t3)'] ^33, 

(10) 



^ [P02 + (Ad - l)Pu] (P22 + P33) + ^ [1 + 2(2rf - 1)(P22 + P33) 



XP33(P12 + P23) - ^ [P02 + P03 + Pl3 + {2d - l)(P2t2 + P3t3)'] P. 
1 



22 



2d 



[1 + 2(2d - 1)(P22 + P33)] P22(Pl2 + P. 



23 J 



^ [Pl3 + (Ad - 1)P23] (P22 + P33) - ^ [P02 + P03 + P13] P33 

+ ^^(At2 + At3)'P33 - ^ [1 + 2(2rf - 1)(P22 + P33)] P33 
X(P12 + P23). 



(12) 



6 



0.25 




).9 0.95 1 



FIG. 2: Stationary order parameter ps versus density in the one-site approximation, for 
d = 1 and d = 2. 



Solution of the above equations is performed numerically. We note that in light 
of the constraints expressed in Eqs. ([1]) and ([2]), we have only five independent 
equations. We also take advantage of the following symmetry: in the mean-field 
approximation, if P02 = -P03 initially, then this equality continues to hold throughout 
the evolution. A similar relation holds between Pi2{t) and Pisit). We therefore 
obtain a set of three independent differential equations for ^22(^)5 P23{t) e P33(t), 
which are readily integrated using a fourth-order Runge-Kutta scheme 23j. We 
define the order parameter as the fraction of active sites. 



(13) 



Numerical integration reveals that in one and two dimensions, p{t) — )■ as t — )■ 00 
for particle densities C ^ 1) while for higher densities it attains a nonzero stationary 
value Ps (see Fig. [2]), which grows continuously with C ~ 1- Thus in the one-site 
approximation, the model exhibits a continuous phase transition between an active 
and an absorbing state. Such a continuous absorbing-state phase transition is famil- 



iar from studies of the contact process 



2J], and of conserved stochastic sandpiles. 



among other models. We verify that for C 7^ p(^) approaches its stationary value 
exponentially: \p{t) — Ps\ oc exp(— t/r), where the relaxation time r depends on 
C, and diverges as C — > Cc = 1, following r ~ 1 / | C — Ccl, as is typical for mean- 
field analysis of absorbing-state phase transitions |25|. The one-site approximation 
yields the critical exponent /3, defined via Ps ~ (C ~ Cc)'^, (for C > Cc), as = 1 
for both d = 1 and d = 2. This value is expected for mean- field analysis of contin- 
uous absorbing-state phase transitions in models that do not possess up-down (or 
particle-hole) symmetry 25j . The reason is that in the absence of such a symmetry, 
all powers of the order parameter are allowed in the mean-field equations of motion, 
so that near the critical point, the terms proportional to p and dominate (that 
is, dp/dt ~ ap — bp"^, with a oc ( — (c and b > 0), and the stationary value of p is 
proportional to C — Cc- 



IV. SIMULATION 



We simulate the restricted sandpile model described above in one dimension us- 
ing periodic boundaries, on rings of L = 500, 1000, 1500, and 2000 sites. The initial 
configuration is defined by assigning limiting heights = 2 or 3 with equal probabil- 
ities, independently, to each site, and then distributing randomly particles among 
the L sites, avoiding occupancies that exceed the maximum height. The resulting 
initial distribution is statistically homogeneous; the occupations of different sites 
are essentially independent. Once all particles have been inserted, the stochastic 
dynamics, which, as noted conserves particles, begins. For each system size L, we 
study an interval of densities ( = N/L. In all cases, we use A^^ = 2000 independent 
realizations of the process. The maximum time is = 5 x 10^ units for L = 1000 
and t„^ = 3x 10^ for L = 2000. 

To determine the critical behavior of the one-dimensional version of sandpile 
defined above, we study the time-dependent density of active sites p{t) as well as 
their survival probability P{t). Figure |3] shows the typical simulation behavior for 
p(t) and P(t). We see that p(t) possesses a transient part before reach a well defined 
stationary value ps{C,L), while the survival probability P{t) oc exp(— t/r(C, L)) has 
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FIG. 3: Simulation: order parameter (upper panel) and survival probability (lower panel) 
versus time for L = 1000 and particle density Q = 1.64. 

an exponential decay. Discarting the initial transient portion of the data the survival 
time r(C, L) is estimated by the slope of the curve. 

In simulations, the particle density C, cannot be varied continuously; for each 
system size L it can only be changed in increments of 1/L. To have acess to intervals 
of particle density smaller than 1/L, we follow a method employed in the study of 
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conserved sandpile models 18| and pair contact process 26|]. Initially we determine 
the stationary average of p for a series of discrete values of the particle density, 
as shown in Fig. HI Since it is reasonable to suppose that the resulting points fall 
on a smooth curve (as is indeed confirmed by the data), we then use quadratic 
interpolation to estimate p at particle densities that are not accessible for the sizes 
studied here. 
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FIG. 4: Simulation: order parameter versus particle density C for sizes as indicated. 



By means the analysis of these data, we obtain the order parameter and the mean 
survival time as functions of system size for diverse values of the particle density, as 
shown in Fig. O At an absorbing-state phase transition, the critical point of a phase 
transition is determined by seeking a power-law dependence of the order parameter 
ps and the survival time r on the system size L. These two parameters are governed 
by 

PsiCL) = L-^/'^-C(L^/^-A), (14) 
t{C,L) = L'^ii/'^^R(L^/''^A), (15) 

where A = — is the distance from criticality and C and R are finite-size scaling 
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relations j27|. At the critical point (A = 0), we expect ps{Q, L) ~ L-l^/^i- and 
r(Cc, L) ~ U\\l^^ . With this in mind, we can estimate Cc and from the curve of 
Ps(Cc, -^^) that best approximates a straight line when plotted versus L on log scales. 
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FIG. 5: Stationary order parameter versus system size L for particle densities C a-s 
indicated. 
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This analysis yields Cc = 1.6400(2) and /3/z/_l = 0.227(5), where the figures in 
parentheses denote the uncertainty in the last significant figure. Analyzing the 
data for the lifetime r in the same manner, we obtain z = = 1.44(3). (The 

uncertainties are related to two contributions: one due to the uncertainty of the 
fit, the other due to the uncertainty in the values of ps and r for each size L.) We 
estimate the critical exponent by plotting L^^'^^ps{(, L) versus L^l^^l\ for various 
system sizes, seeking the value of which yields the best data collapse. Figure [6] 
shows that a good collapse is obtained using 1/z/_l = 0.704(5). The critical exponent 
/3 that relates the order parameter pg with A through the relation ps = is then 
easily determined as /3 = 0.322(5). 

Table [T] compares our estimates for critical exponents with those obtained in 
studies of other one-dimensional models in the CDF universality class. Despite 
apparent differences, it is important to recall that previous studies have revealed 
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FIG. 6: Scaling plot for the density of active sites. 
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TABLE I: Critical exponents for one-dimensional models in the CDP universality class compared 
with estimates from present work. 
Sleepy Random Walkers 



Restricted Manna model 



3C 



28|; CDP Field theory 



29 



Model 



Rest. Manna^ 0.213(6) 1.55(3) 0.29(1) 
CDP - FT"^ 0.214(8) 1.47(4) 0.28(2) 
SRW^ 0.212(6) 1.50(4) 0.290(4) 

Present work 0.227(5) 1.44(3) 0.322(5) 



that simulations of large systems (20 000 sites or larger) are needed to obtain reliable 
values of critical exponents for this class 28|, l30 |. 



V. CONCLUSIONS 



We study a height-restricted fixed-density version of the Oslo sandpile in one 
dimension. At each site, the limiting height Zc may be either 2 or 3. The model is 
found to exhibit a continuous phase transition between and active and an absorbing 
state at a critical value of the particle density, Cc- The one-site mean- field approx- 
imation predicts Cc = 1, whereas simulations yield Cc = 1.6400(2). The small sizes 
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analyzed here limit the reliability of our estimates for the critical exponents. Given 
the comparison with the literature values shown in Table 11, a tentative conclusion 
that the restricted Oslo model belongs to the conserved directed percolation class 
seems reasonable. More definitive conclusions will require studies of larger systems. 
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